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ABSTRACT 

We present numerical radiation-hydrodynamic simulations of the evolution of HII regions formed in an 
inhomogeneous medium resulting from turbulence simulations. We find that the filamentary structure of the 
underlying density distribution produces a highly irregular shape for the ionized region, in which the ionization 
front escapes to large distances in some directions within 80,000 years. In other directions, on the other hand, 
neutral gas in the form of dense globules persists within 1 parsec of the central star for the full duration of our 
simulation (400,000 years). Divergent photoablation flows from these globules maintain a root-mean-squared 
velocity in the ionized gas that is close to the ionized sound speed. Simulated images in optical emission lines 
show morphologies that are in strikingly detailed agreement with those observed in real HII regions. 

Subject headings: HII regions — ISM: clouds — ISM: kinematics and dynamics — stars:formation — turbu¬ 
lence 


1. INTRODUCTION 

The currently accepted view of the interstellar medium 
(ISM) is that its density an d velocity distributions are s haped 
by turbulence (see , e.g. Vazquez-Semadeni et al. 2000 


jtlmegreen & Scalo 2004). In the molecular clouds that will 


constitute the sites of star formation, the ISM turbulence is 
supers onic, with tu rbulent Mach numbers as large as 10 or 
more (|^arson| 1981). Such high Mach numbers indicate that 
the medium is highly compressible. The efficiency of the star 
formation process is known to be small from observations, 
at mo st 10%-30% for cluster-forming cores (Lada & Lada 
2003). In numerical simulations of supersonic turbulent me¬ 
dia, large-amplitude density fluctuations are generated, which 
may become locally gravitationally unstable and collapse to 
form stars. 

Once a massive star forms, it will photoionize its surround¬ 
ings, forming an H n region. Observed H II regions display a 
wide variety of morphologies and sizes and are generally ir¬ 
regular. While still embedded in their natal molecular clouds, 
HII regions remain small and are classified as ultracompact 
(lin earsize<0.1 pc) or compact (0.1-1 pc) ( WV>od & Church- 
well 1989| ). Extended HII regions (size > 1 pc) are thought to 
correspond to more evolved states. The formation of HII re¬ 
gions and the propagation of ionization fronts in uniform me¬ 
dia has been well studied (Stromgren 1939 |<alin| |l 954| ). The 
expansion of HII regions in power-law and plane-parallel den¬ 
sity distributions, and the breakout of ionized gas from dense 
clouds into a low density medium has also been in vestigated 
( Tenorio-Tagle|l979|: |Bedijn & Tenorio-Tagle|l9~8l1; Franco e t 


al. |1989[ |1990| ; |Henney et alj|2005a| : Arthur & Hoare||2006| ). 

These models all assume a smooth density variation and do 
not give rise to the irregula r complex shapes exhibited by real 
HIIregions (Henney 2006). 
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It has been suggested that the irregular shapes of observed 
Hu regions, including the “fingers” and “elephant tru nks” 
seen in regions such as M16 (see, e.g., ^-tester et~ak 1996), are 
due to instabilities in the ionization front. When cooling is in¬ 
cluded in radiation-hydrodynamical simulations of HII region 
expansffiiLa^tron^mstabilit^j^redicte^anal^ticall^h^Giu- 


liani ( 1979| ), le ads to vigorous fragmentation o f the expanding 
massive shell (Tarcla-Segura & Franco |l 996). When small- 
scale den s ity flu ctuations are present in the ambient medium, 
Williams (1999) shows that a shadowing instability in the R- 
type ionization front phase can lead to a non-linear hydrody¬ 
namic instability and the formation of clumps in the D-type 
phase. The clumps and filaments formed by such instabili¬ 
ties tend to be radial in nature, and thus do not account for 
observed structures such as the Orion bar, which has a linear 
appearance, perpendicular to the direction to the ionizing star. 

It is entirely possible that the irregular appearance of ob¬ 
served HII regions is due not to instabilities in the ionization 
fron t but to underlying structure in the ambient medium. L i 


et al. ( 2004 ) calculated the propagation of the R-type ioniza- 
tion front produced by a massive star in a precomputed three- 
dimensional compressible turbulent density field. Their re¬ 
sults show that the initial density structure is important for 
the resulting structure of the H II region. However, these cal¬ 
culations do not include hydrodynamics and therefore can¬ 
not follow the subsequent dynamical expansion of the ion¬ 
ized gas into the surrounding tu rbulent med i um af ter the ini¬ 
tial ~ 100 yrs. More recently. Dale et al. (2005) have pre¬ 
sented smooth-particle-hydrodynamic simulations of the de¬ 
velopment of an H II region around a newly-formed star clus¬ 
ter. These simulations do follow the dynamical expansion of 
the ionized gas in a highly clumped medium and its feedback 
on the star formation process, but the initial density field is 
specified in a somewhat ad hoc manner. 

In this article we present results from the first fully coupled, 
radiation-hydrodynamic simulations of the formation and ex¬ 
pansion of HII regions in a turbulent medium. As our starting 
point, we take a precomputed three- dimensional compressible 


t urbul ent density field calculated by [Vazquez-Semadeni et al. 
(2005). In § pf we describe our numerical method and the 


initial conditions. Our results are presented in § |] in a form 
directly comparable to observations. In § ^ we summarize and 
conclude our findings. 
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2. SIMULATIONS 
2.1. Numerical Method 

The hydrodynamic s is calculated using th e non- relativistic 
scheme described in ^uldcrink & Mell ema (1995 ), with the 
addition of a Local Oscillation Filter (Sutherland et al. 2003) 


to suppress numerical odd-even decoupling behind radiatively 
cooling shock waves. For the radiative transfer, we implement 
t he ne w, efficient C 2 -Ray method, described by Mellema et al 


(|200fJ), which is explicitly photon conserving. The C T -Ray al 
gorithm uses an analytical relaxation solution for the ioniza¬ 
tion rate equations, which allows for time steps much larger 
than the characteristic ionization timescales and the timescale 
for the ionization front to cross a numerical grid cell. This 
efficiency allowed us, in another current application, to use 
this code to perform the first large-scale simulations of the 
reionization of the universe, which required solving for the 
radiative feedback of up to tens of thousands of s ources on 
200 3 -400 3 computational grids ( Iliev et al. 2006a ). The ra¬ 
diative transfer is coupled to the hydrodynamics via opera¬ 
tor splitting and tests of the com bined numerical method are 
presented in |Hiev et al. (2006b, B[). 6 These features make the 
current problem of three-dimensional HII region expansion 
in a turbulent medium entirely tractable, and reduce the com¬ 
putational time for the simulations to a few hours for a 128 3 
grid on a single processor. The simulations presented in this 
paper are for a single, fixed grid since the radiative transfer 
algorithm is parallelized only for shared-memory machines. 


2.2. Initial Conditions 

As initial conditions for the density and velocity fields we 
take the three-dimensional numerical simulations of driven 


turbulence presented by Vazquez-Semadeni et al. (2005). 


These turbulence models are isothermal but once included in 
our simulations the gas is subject to h eating due to photoion¬ 
ization, and cooling, calculated as in Raga et al. (1999). We 
use the simulations with root-mean-squared (rms) sonic Mach 
number M s = 10, which gi ve a large dynamic range for the 
initial density fluctuations. Vazquez-Semadeni et al. (2005) 
present both ideal-magnetohydrodynamics (MHD) and non¬ 
magnetic cases but we have chosen a nonmagnetic simula¬ 
tion for the present paper because our hydrodynamics scheme 
does not include MHD. The turbulence simulations are scale 
free and are characterized by three nondimensional numbers: 
Mj = cr/c (the rms sonic Mach number, where cr is the turbu¬ 
lent velocity dispersion and c is the sound speed), J = L/Lj 
(the Jeans number, giving the size of the box in units of the 
Jeans length Lj), and (3 = Fth/Rmag (the ratio of thermal to 
magnetic pressures). In the case we consider, M s = 10, J = 4, 
and (5-oo. The turbulence calculation is time dependent 
and we choose a late time, when several regions are collaps¬ 
ing, as the initial state for the principal run of our radiation- 
hydrodynamics simulation. 

The scaling to physical variable s we choose is that used 


by Vazquez-Semadeni et al 


(2005), that is, isothermal sound 


speed c s = 0.2 km s 1 and mean nucleon number density 


6 Details of the test pro hkms can he found aj_j ittp://www. 
mpa-garching.mpg.de/tsu3A 

7 The original simulation used molecular gas with mean molecular mass 
2.4»7 h and number density 500 cm -3 . Since we do not treat the dissociation 
of molecular hydrogen in our code, we have substituted atomic hydrogen 
with the same mass density and sound speed. Another consequence of this 
change is that our initial temperatures are lower, but this is of no dynamical 
consequence. 


no - 932 cm~ 3 . This gives a spatial size L - 4 pc for the com¬ 
putational grid. We place the ionizing photon source at the 
center of the densest clump, which has density n = 1000«o, 
and the periodic boundary conditions used for the turbulence 
simulation enable us to move this position to the center of the 
computational grid. We found that the densest clump in the 
turbulence simulations has a large absolute velocity, since it 
forms due to velocity fluctuation in the molecular gas. Con¬ 
sequently, it is necessary to shift the velocities of the whole 
grid by an amount equal to the velocity of the densest cell, in 
order that the densest material remains at rest in the frame of 
the star. For the ionizing source, we choose a stellar ioniz¬ 
ing photon rate 5, = 5 x 10 48 s 1 and effective temperature 
Tpff = 37500 K, equivalent to an 07.5 main sequence star 
( Panagia 1973 ). In order to account for the material that goes 
into the formation of the ionizing star, we remove 90% of the 
mass from the central 2x2x2 group of cells, which corre¬ 
sponds to = 27M e . 

Although the initial turbulence simulation has periodic 
boundary conditions, our hydrodynamics calculation has open 
boundaries. This is because in cases where the photoionized 
region expands beyond the grid boundary, periodic bound¬ 
ary conditions would be unphysical. Furthermore, our cal¬ 
culations do not include self-gravity and our treatment of the 
heating and cooling of the cold neutral/molecular gas is only 
approximate and does not include chemistry. 


3. RESULTS 

The evolution of the ionized region is illustrated in Fig¬ 
ure [I] by means of density and temperature slices in different 
yz planes for a sequence of times. At t - 50,000 years, the 
ionized gas is confined to a small volume around the star, as 
can best be appreciated in the temperature slices. One also 
sees high temperature gas behind the shock that precedes the 
ionization front, and which has now detached from the ion¬ 
ization front as it encounters lower density gas. Between the 
shock and the ionization front is a high density swept-up shell 
of neutral gas. By t = 100,000 years, the ionized region has 
expanded considerably, reaching a radius of ~ 1 parsec in 
most directions and even escaping from the grid in one direc¬ 
tion. Some dense neutral gas remains close to the ionizing 
star, however, which forms photoablated globules and fila¬ 
ments. The transonic flows from these can be seen to evacuate 
cavities in the ionized gas, bounded by shocks. At later times, 
the photoablation flows become even more prominent as more 
dense neutral clumps are uncovered, while existing globules 
are eroded and pushed away from the star by the rocket ef¬ 
fect. By t - 400,000 years the ionization front has escaped 
from the grid in most directions but some neutral gas remains 
as close as 1 parsec to the ionizing star. The temperature is 
roughly uniform in the ionized region at = 8000 K, reaching 
slightly higher values (= 10,000 K) close to the ionization 
front and in shocks, and slightly lower values (= 7000 K) due 
to expansion cooling in the photoablation flows. 

Figure [l] (top panel) shows how the ionized fractions of 
the total volume, X vo i, and mass, X mass , of our simulation 
varies with time. Both fractions are very small for the first 
= 50,000 years of evolution, during which the ionization 
front remains trapped inside the dense clump in which the 
star formed. By the end of our simulations, over 70% of the 
volume of our cube has been ionized, but only 15% of the 
mass. However, this latter figure is only a lower limit, since it 
only accounts for the ionized gas that remains on the grid (see 
below). 













































HII regions in turbulent clouds 


3 





Eh 



-1 0 +1 

x offset from star, parsec 


Fig. 1.— Gas density and temperature slices through our simulation for a sequence of evolutionary times, increasing from top to bottom, as indicated on the 
left axis. Density is shown on a negative logarithmic scale between 10 (white) and 10 5 cm -3 (black). Temperature is shown on a positive linear scale between 0 
(black) and 10 4 K (white). For each time, a sequence of 8 slices in the yz plane are shown at x intervals of 0.25 parsec, each with size 4 by 4 parsecs (x offsets from 
the ionizing star position are indicated on the bottom axis). This figure is also available as a color mpeg animation in the electronic edition of the Astrophysical 
Journal. In the animation, the density of cold neutral gas is shown in positive grayscale, warm neutral gas in blue/green, and ionized gas in red/orange for three 
perpendicular midplane slices. 
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Fig. 2. — Top: volume-weighted (solid line) and mass-weighted (dashed 
line) ionized fraction of the computational grid against time. Bottom: evolu¬ 
tion with time of the mean radius (solid line) of the ionized region, compared 
with the minimum (dashed line) and maximum (dotted line) radius of the ion¬ 
ization front at each time. Also shown is the analytic solution (thick gray line) 
for the evolution of the mean radius of a Stromgren sphere in a homogeneous 
medium with the same mean density as our simulation. 


The lower panel of Figure 2] shows the evolution of the 
mean radius of the ionization front, averaged over all direc¬ 
tions from the ionizing source, (f?; on ) = (32^ _ vo lVg r id/47^) 1/,3 , to¬ 
gether with the minimum and maximum ionization front radii 
at each time, R m j n and R max , respectively. The ratio R m n*/ Rmm 
increases with time, indicating an increasingly irregular shape 
for the ionized region, until at / ^ 80,000 years the ion¬ 
ization front suddenly escapes from our computational box 
in some directions, at which point /? max becomes limited by 
the size of the grid. R mm remains small at < 1 parsec dur¬ 
ing the entire simulation due to the survival of dense neu¬ 
tral clumps as discussed above. For comparison, we also 
show the evolution of a D-type ionization front in a uni¬ 
form density medium with the same mean density as our 
simulations, which is given by the expression (e.g., [Spitzei 


19681) R = R 0 (l + jCit/Ro) A/1 , where Ci is the ionized sound 


speed and Rq = (S= 0.532 parsec is the initial 
Stromgren radius (the initial R-type propagation of the front 
lasts for only - 100 years for our parameters). At early times, 
our radii are far below this prediction since the 1000 times 
overdensity at the position of the ionizing star gives an initial 
Stromgren radius that is 100 times smaller. By 150,000 years, 
the mean radius has caught up with the uniform case, at which 
point its expansion appears to slow, although this is probably 
simply an effect of the finite size of our grid. The minimum 
radius is always well below the uniform case and this is in¬ 
sensitive to the size of the grid. 




Fig. 3. — Top: Mean density, (n) against time. Bottom: Clumping, C = 
( n 2 )/(n ) 2 , against time. Solid line — ionized gas; dashed line — neutral gas; 
dotted line in bottom panel — reciprocal of filling factor, s = (/7 2 ) 3 /(« 3 ) 2 , of 
ionized gas. 


Figure £?] shows the evolution of the density in our simula¬ 
tion. The mean density (top panel) in the ionized gas shows 
a dramatic decline as the Hu region expands, whereas the 
mean density of the neutral gas shows a slight increase af¬ 
ter 100,000 years due to the effects of shocks driven by the 
expanding ionized region. The slight decline in the mean to¬ 
tal density towards the end of the simulation indicates that 
we are starting to lose mass from our grid. Although we do 
not track the ionization state of this gas, we expect it to be 
predominantly ionized, giving an upper limit to the fraction 
of the initial simulation mass that has been ionized by the 
end of the run of - 50%, as opposed to the lower limit of 
15% that was derived above. The density clumping of the 
neutral gas (bottom panel) is initially very high due to the 
gravitational collapse of dense regions in the turbulence sim¬ 
ulations but quickly falls with time as the densest region be¬ 
comes photoionized and expands. The clumping of the ion¬ 
ized gas is much less pronounced, with (n 2 )/(n) 2 not exceed¬ 
ing 2. However, higher order statistics indicate significant in¬ 
homogeneities in the ionized gas too, with a filling factor of 
e ~ 0.02 during much of the evolution. This is primarily 
due to photoablation flows from globules, which have ionized 
density contrasts of 3-10 with respect to the mean, and partly 
due to weak shocks, which have typical density contrasts of 
2 ^ 1 . 

Figure |4 shows the mean gas velocities in the simulations. 
The one-dimensional rms velocity over the volume, V , of our 
simulation is defined as 

<y2>1/2 = (fe |V|2 dV 1 3 fe dV ) ’ 
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Fig. 4.— Average gas velocities against time. Solid line — one-dimensional 
rms velocity of ionized gas; dashed line — one-dimensional rms velocity 
of neutral gas; dotted line — mean expansion velocity of ionized gas; dot- 
dashed line — mean expansion velocity of neutral gas; thick gray line — 
mean expansion velocity of ionized gas in a Stromgren sphere evolving in a 
uniform medium with the same mean density as our simulation. 


and the mean expansion velocity as 


< Vr) 


I^ idV l f‘ 


fdV, 


where v is the vector gas velocity at a point, r is the unit vector 
in the radial direction away from the ionizing star, and £ is the 
ionized or neutral fraction, as appropriate. The rms velocity 
of the ionized gas increases with time for the first 50,000 years 
until it is close to the ionized sound speed, from which point 
on it remains roughly constant. The rms velocity of the neu¬ 
tral gas remains at its initial value of ~ 2 km s -1 for the first 
100,000 years, after which it is increased to ^ 5 km s 1 by 
the action of shocks driven by the expanding ionized gas. 
The mean neutral expansion velocity is a significant fraction 
of the rms velocity, indicating that the neutral gas motions 
are predominantly away from the ionizing star. The same is 
true for initial expansion of the ionized gas, but after 100,000 
years the ionized expansion velocity gradually declines to 
only 2 km s' 1 , while the rms velocity remains at 8 km s' 1 . 
This is due to the effect of the photoablation flows from dense 
globules, which flow back towards the star at speeds of up to 
20 km s' 1 and partially cancel out the general expansion of 
the ionized gas. The mean expansion velocity of the ionized 


gas in a uniform medium is equal to | dR/dt and is always less 
than half of the ionized sound speed, which is much lower 
than is seen in our simulation. 

Figure || shows 3-color optical emission line images of 
our simulations at three different evolutionary times, calcu¬ 
lated using 2-level at om approxima ti ons fo r the line emis- 
sivity as described in |Henney et al. ( 2005t ). The ion frac¬ 
tions of N and O were approximated as functions of the H 
ioni zation fraction, which were calibrated using Cloudy (F er- 


land 200C ). This is a reasonable approximation since depar¬ 
t ures from equilibrium ionization should be similar for all ions 
(Henney et al. 2005b). The images include the effects of local 


absorption by dust in the neutral gas, w ith an assumed cross 
section of 5 x 10' 22 cm 2 per H nucleon (Baldwin et al. 1991) 
but do not include scattering. They are designed to be di¬ 
rectly comparable with observed narrow-band emission line 
observations of real nebulae. In panel a, although the ion¬ 
ization front has already broken out in some directions, this 
is only where the density is very low, resulting in weak emis¬ 
sion. Bright cusps indicate the ionization fronts at the heads of 
photoablated globules that protrude into the interior of the ion¬ 
ized region. On a larger scale, the outer border of the emission 
seems to consist of several roughly straight segments. This is 
due to the formation of cylindrical ionization fronts around 
filaments in the neutral density distribution. Dark lanes are 
due to dust absorption in dense filaments that lie outside the 
ionized region. As the simulation progresses (panels b and c), 
one sees a greater number of photoablated globules and at the 
same time, the overlying extinction becomes less. Shocks due 
to the collisions between photoablation flows produce bright¬ 
ness variations in the fully ionized gas. 

4. DISCUSSION AND CONCLUSIONS 

HII regions around single stars are found with a variety of 
sizes, from ultracompact (< 0.01 parsec) to extended (> 1 par¬ 
sec) and it is probable that these represent an evolutionary se¬ 
quence (e.g., jGarcfa-Segura & Franco 1996 ). Our simulation 
passes through all of these phases during its evolution but, be¬ 
cause of the physical scaling we have chosen, the initial com¬ 
pact phases are very poorly resolved. We therefore compare 
our results with the observed structure and dynamics of more 
evolved, extended H II regions. The shapes of such regions 
tend to be very irregular on a large scale and frequently show 
bright-rimmed structures on smaller scales ( Pottasch 1956 ), 
together with filamentary overlying extinction features. Many 
Galactic H II regions show strong similarities in appearance 
with our simulation (see especially Fig. 5|), with perhaps the 
most striking rese mblance occurin g in the case of M20, the 
Trifid Nebula (see Rho et al. 2005 ). In the case of many other 
regions, such as the Orion Nebula (M42), and the Eagle Neb¬ 
ula (Ml6), the details of the observed emission structures are 
very similar to those seen in our simulations but the large- 
scale distributions are somewhat different. This is presum¬ 
ably due to density gradients on a scale equal to or larger than 
that of the H II region itself, which give rise to champagne 
flow s that are ionization-bounded on one side (Boden heimer 
et al. 19791; [Tenorio-Tagle||1979|; Henney et al.||2005~al; Arthur 


& Hoare [2006[ ). These are not present on our simulation, in 
which the density is roughly homogeneous on scales > 1 par¬ 
sec, although the effects of such l arge-scale g r adient s can be 
clearly seen in the simulations of pale et ah (2005). On an 
even larger scale, very similar emission structures are seen in 
giant H II regions, of which NGC 3603 in our Galaxy ( M()flat 
et al. |2003) and 30 D oradus in the Large Magellanic Cloud 
( {Scowen et al. 1998) are the best-studied examples. How¬ 
ever, such regions are ionized by multiple stars and therefore 
cannot be directly compared with our si mulations, although 


simila r processes are expected to occur (I’enorio-Tagle et al 

2006b. 


There are hints in our simulation of a qualitative change in 
the morphology of the region as it evolves. At times soon after 
the breakout of the ionization front from the dense core, one 
finds photoablating structures that are predominantly cylin¬ 
drical in form, oriented side-on to the flux of ionizing pho¬ 
tons. These give rise to bar-like features in the emission maps, 
which are similar in appearance to the Bright Bar and other 
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Fig. 5.— Synthetic narrow-band optical emission line images in the light of [N II] 6584 A (red), Hir 6563 A (green), and O HI 5007 A (blue) for our simulation 
at evolutionary times of (left panel) 100,000 years, (center panel) 250,000 years, and (right panel) 400,000 years. In each case, projections are along the six 
principal axes of the simulation—top row: +x, -x\ middle row: +y, -if. bottom row: +z, — z. The scales are linear in surface brightness, with maximum values 
of 0.005 ([N II] 6584 A), 0.028 (Ho- 6563 A), and 0.018 (O III 5007 A), all in units of erg s' 1 crrT 2 sr~' and assuming gas-phase abundances of O/H = 3 X 10 -4 
and N/H = 5 X 10 -5 . 


simila r features found in the Orion Nebula (O ’Dell & Yusef- 
Zadeh 2000 ). At later times, the photoablating structures be¬ 
come more finger-like, with a spherical ablation flow from the 
bright tip of the finger, which points towards the ionizing star. 
An additional feature of the late-time evolution of our simula¬ 
tion is the formation of recombination fronts, which temporar¬ 
ily reverse the expansion of the ionized region in certain di¬ 
rections. These form due to collisions between ionized flows, 
forming dense shells that in certain circumstances are capable 
of trapping the ionization front. Although these recombina¬ 
tion fronts also show finger-like irregularities, they are easily 
distinguished from the photoablating fingers because the ion¬ 
ized density does not show a maximum at the ionization front. 
As a result, they give soft, diffuse edges to the H II region in 
the emission maps, as opposed to the sharp edges associated 
with the photoablation flows. 

The rms velocity of the ionized gas is remarkably constant 
at ^ 8 km s 1 after the initial breakout of the ionization front 
from the dense core, which occurs around t - 50,000 years. 
For the first ^ 150,000 years, this is driven largely by the 
physical expansion of the ionized core in a mode similar to 
the full^_densit^bounded_cham£agne^ows_studied^^jShu 
et al. ( |2002| ). At later times, it is the divergent photoablation 
flows from numerous neutral globules that maintain the rms 
velocity. The rms velocity of the neutral gas is also strongly 
affected by the ionized region, increasing to ^ 5 km s _1 af¬ 
ter 150,000 years, which is 2-3 times the initial turbulent 
value. Optical spectroscopy of Hu regions frequently shows 
non-thermal linewidths of ~ 10 kms _1 , which are roughly 
constan t betw e en partially ionized a nd fully-ionized species 
p’Dell| (2001); O’Dell et al. (2003). This is difficult to ex¬ 
plain in the context of HII regions in smooth density gradients 


Our simulations are only a first step towards fully self- 
consistent models of H n region evolution in realistic density 
distributions. Important physical processes that we have ne¬ 
glected include the stellar wind from the ionizing star and the 
radiative acceleration of dust grains, which are coupled to the 
gas via Coulomb collisions. Both these would tend to evacu¬ 
ate a central cavity in the photoionized gas. Even in smooth 
density distributions the effect of stellar winds on the dynam¬ 
ics of th e ionized gas is found to be rather complex ( Arthur 
& Hoare 2006| ). The effects of magnetic fields on the ioniza- 
tion fronts and on the dynamics of the neutral gas (Redman 


( Henney et al. 2005a ) but would be a natural consequence of 
our simulation. 


et al. |1998| ; |Williams et al.| |2000| ; |Williams & Dyson 2001) 
is another area that needs to be explored, since it is possi¬ 
ble that MHD effect s partially dictate the shape of the pho- 
toablated struct ures (Carlqvist et al. 2003; Ryutov et al. 2005; 
William s|p006| ). A more realistic treatment of the diffuse ion¬ 
izing field may also affect the details of the shapes of the pho- 
toablated globules but should not significantly change the dy¬ 
namics. The evolution of the region on timescales longer than 
the 400,000 year duration of our simulation is also of interest, 
but would require the consideration of a larger computational 
domain. 

In summary, we have carried out the first investigation of 
the hydrodynamic evolution of an H II region in a realistic tur¬ 
bulent medium. We have shown that the clumpy and filamen¬ 
tary structure of the underlying density distribution leads to an 
extremely irregular shape for the ionization front. The den¬ 
sity structure in the ionized region is less clumped than that 
of the neutral gas, but nevertheless shows complex dynam¬ 
ical structures due to the mutual interactions between pho¬ 
toablation flows from neutral globules and filaments, some of 
which survive at close distances to the ionizing star (< 1 par¬ 
sec) throughout the evolution. The rms velocity of the ion¬ 
ized gas is approximately equal to the ionized sound speed 
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(~ 10 km s 1 ) during almost the entire evolution, even at late 
times when the net radial expansion is very low (< 2 km s -1 ). 
The calculated emission-line morphologies of our simulation 
show striking similarities to the appearance of real Hll re¬ 
gions. 
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